
module IO

implicit none

contains

 Function numcat(S, num)
   character(LEN=*) S
   character(LEN=120) numcat, numstr
   integer num

   write (numstr, *) num
   numcat = trim(S) // trim(adjustl(numstr))
   !OK, so can probably do with with a format statement too... 
 end Function numcat


 Function sortedTauIntensity(ntau, values)
  integer,intent(in)      :: ntau
  real   ,intent(in)      :: values(:)
  real                    :: tauIntensity(2*ntau)
  real                    :: sortedTauIntensity(2*ntau)
  real, dimension(:)      :: tau(1:ntau), indexMax(ntau)
  real                    :: tauMax 
  integer                 :: ii, jj
  tauIntensity(1:2*ntau-1) = values(1:2*ntau-1)
  tauIntensity(2*ntau) = 1.00 
  do ii = 1, ntau-1
        tauIntensity(2*ntau) = tauIntensity(2*ntau) - values(2*ii)
        tau(ii) = values( 2*ii-1 ) !int(100*values( 2*ii-1 ) )
  enddo
  tau(ntau) = values( 2*ntau-1 ) !int(100*values( 2*ntau-1 ) )
  !print*, "---------------------"
  !print*, tauIntensity
  !print*, "tau=", tau
  do ii = 2*ntau-1, 1, -2 
        !do jj = 1 , ntau
        !tauMax = 
        indexMax = MAXLOC(tau)
        !write(*,*) "ii=", ii, ",  indexMax=",indexMax(1) , ", tau(indexMax)=", tau(indexMax(1) )
        sortedTauIntensity( ii:ii+1 ) = tauIntensity( 2*indexMax(1)-1:2*indexMax(1) ) 
        tau(indexMax(1) ) = 0.0
  end do
  !print*, "tau=", tau
  !print*, sortedTauIntensity
  !print*, "---------------------"
 end Function 
 
end module IO


program main 
 !use  myprecision
 use  IO
 implicit none
 character(LEN =50) :: fmt
 integer,parameter   :: ntau=5
 real                :: values(1:2+ 2*ntau-1) = (/0.99,0.99, 0.55,0.1,  0.44,0.2,  0.22,0.3,  0.33,0.2,  0.11 /)
 real                :: main_tmp(2*ntau)

 fmt = trim(numcat('(2E16.5,',ntau*2-1))//"E16.7,"//'E16.5)'
 print*, "(fmt: ", fmt,")"
   
 print trim(numcat('(A,',size(values) ))//'E16.7,A)', "(originalTauIntensity: ", values,")"
 main_tmp = sortedTauIntensity(ntau, values(3:2*ntau+2) )
 print trim(numcat('(A,',ntau*2+2))//'E16.7,A)',      "(sortedTauIntensity  : ", values(1:2), sortedTauIntensity(ntau, values(3:2*ntau+2) ), ")"

end program 

